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Sine-based method for a high-precision solution in the direct space of wave equations 

with periodic boundary conditions 

Paolo Marconcini/ '|j Demetrio Logoteta/ and Massimo Macucci-'^ 
^ Dipartimento di Ingegnena dell'Informazione, Universitd di Pisa, Via Caruso 16, 1-56122 Pisa, Italy. 

Direct-space finite-difference approaclies usually require, to achieve a similar accuracy, the solution 
of larger algebraic problems with respect to those resulting from a formulation based on the Fourier 
transform in the reciprocal space. This is the result of errors that are introduced in the treatment 
of the derivatives with the application of direct-space discretization formulas. Here we propose an 
approach, based on a particular set of basis functions, that allows us to obtain an exact treatment 
of the derivatives in the direct space and that is equivalent to the solution in the reciprocal space. 
We apply this method to the numerical solution of the Dirac equation in an armchair graphene 
nanoribbon with a potential varying only in the transverse direction. 

I. INTRODUCTION 



With the progressive downscaling of electron devices, whose size has reached well into the nanoscale, the quantum 

simulation of transport has become increasingly important. One of the key points in the solution of the differential 

C^ ' equations describing charge transport is the choice of the solution domain. Indeed, in the case of periodic boundary 

rS^ ' conditions, the solution in the Fourier domain represents a valid alternative to a study in the direct space. In a 

Jy^ . Fourier analysis the functions appearing in the equation arc replaced by their Fourier coefBcients and the unknowns 

QJ ' are represented by the Fourier coefficients of the wave function (instead of by the values at the nodes of a discretization 

grid, as in direct space approaches). 

The main advantage of this technique is the correct treatment of the derivatives. In a reciprocal space analysis 
C^ I each of the complex exponential functions appearing in the Fourier expansions can be derived exactly, without 
any approximation. On the contrary, in a standard finite difference solution (in the direct space), the derivatives 
are substituted by finite difference approximations involving a certain number of samples of the function. This 
■ O . approximation introduces a distortion in the dispersion relation, which severely limits the accuracy of high-order 
solutions. The discretization error decreases if finer grids are used, at the expense of an increase of the computational 
effort and of a consequent reduction of the code performance. Moreover, in some cases, such as the solution of the 
Dirac equation for massless fermions, the discretization of the equation on a direct space grid gives rise to problems 
like the appearance of spurious solutions or fermion doubling, unless proper, non-standard discretization techniques 
are applied^"— , while these issues do not appear using a Fourier solution technique. 

Here wc illustrate how the use of a particular set of periodic basis functions makes it possible to obtain a correct 
description of the derivatives in the direct space and a solution approach exactly equivalent to that in the reciprocal 
psj ' space. 

■^ , After discussing how each term of the equation should be treated, we apply this method to the solution of the Dirac 

equation in an armchair graphene nanoribbon with a potential varying only in the transverse direction. Graphene 
(^ I represents a novel material that, since its isolation in 2004^, has been the focus of a vast theoretical and experimental 
(^ • research activity^i^. Among the many applications that have been proposed for this material, there is the possibility 
to use graphene nanoribbons for the implementation of nanoelcctronic devices. The solution of the envelope function 
equation (which in the case of graphene corresponds to the Dirac equation^) within a slice with a potential varying 
only in the transverse direction can be exploited to perform a transport analysis for an armchair nanoribbon with a 
generic potential, if wc ideally subdivide the ribbon into a series of slices in each of which the potential is approximately 
constant in the longitudinal direction and then we apply a recursive technique to obtain the transmission in the overall 
device^. 

II. MATHEMATICAL METHOD 

Let us consider a wave equation defined on a domain [0, L] and with periodic boundary conditions. 

For the sampling theorem, if a function z{x) is assumed to be band limited with band B, it can be perfectly 
reconstructed from its samples taken with a sampling interval A, if i? < 1/(2 A) (Nyquist criterion), according to the 
following equation^: 

<^)= Y. zinA)smcf^^^) , (1) 



where the sine function is defined as sinc(a;) = sin(7rx)/(7rx). 

Moreover, if z(x) is periodic with period L and we take N samples within a period (and thus L ~ N A), we have 
that 

= X:^(^A) g sinc( "-(' + ^^)^ )^X:'^(^A)g^^A(x), (2) 

where we have exploited the periodicity of z[x) and we have defined the function (periodic with period L = N A) 
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If we define the scalar product between two functions "01(2;) ^nd 'ip2{x) as 
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we can prove that the set of functions gt^ix) are orthonormal, i.e. {gj,Aix)\g£A{x)) = Sj^e. 

In our method we express each of the terms appearing in the wave equation as a linear combination of the functions 
gt.,A{x), with coefficients given (in their turn) by linear combinations of the samples of the (unknown) wave function 
p(x) at the points n A (with n = 0, . . . , A^ — 1). In order to easily arrive at closed-form analytical expressions, we 
assume A'' to be an odd number, i.e. iV = 2D + 1 (with D an integer). We rewrite the wave equation in this form, 
project it onto the functions Qj^ix) (with j = 0, . . . , A^ — 1), and exploit the orthonormality of the g functions: we 
thus recast it into a system of N linear equations, that in the following section we will prove to be exactly equivalent 
to a Fourier analysis with a frequency cut-off. 

We will consider the following terms: a term proportional to the unknown wave function p(x), a term proportional 
to a derivative of the wave function, and a term proportional to the product between the wave function p{x) and a 
function representing the behavior of the potential energy inside the domain. Other possible terms appearing inside 
the wave equation will undergo a similar treatment. 

III. TREATMENT OF THE WAVE FUNCTION AND OF ITS DERIVATIVES 



We assume the unknown periodic wave function p[x) to be a band limited function with a maximum band B = D/L, 
such that, sampling it with a step A, the Nyquist criterion B < 1/(2 A) is satisfied. 

Therefore the term proportional to the unknown wave function is easily rewritten in the desired form by sampling 
it with a step A: 
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p{x)= Y,p{lA)ge^A{x). 



(5) 
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Concerning the term proportional to the derivative of the wave function p{x), if p{x) is band limited with a band 
equal to B = D / L also its derivatives have the same property, and thus we can apply the sampling theorem to them, 
too, obtaining (for the generic s-th derivative): 
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We now have to express the values of the derivatives in the points /i A as a function of the samples of p{x). This can 
be easily done by deriving Eq. ([5]). We have that 
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where we have defined as a*^ the coefficients with which we have to combine the samples of the unknown wave 
function p{x) in order to find the samples of its s-th derivative. 

In particular, let us develop the calculation for s = 1 (first derivative) and for s ~ 2 (second derivative). 

For s = 1 we have that 
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and thus 
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If we develop the calculation assuming N odd, we have that 
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Analogously, for s = 2 we have that 
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and thus 
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Assuming TV odd, we derive that 
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Notably, these expressions coincide with those of the so-called SLAC derivative techniquoi^, that have been previ- 
ously obtaineclii switching to the reciprocal space, operating in that domain and then transforming back (analogously 
to what we do in Section |Vl where we compare our method with the Fourier one). 



IV. TREATMENT OF THE PRODUCT TERM 



The last terms that we have to manage are those proportional to the product between the potential energy function 
f{x) and the wave function p{x). 



We can assume, as a first approximation, that also this term is band limited, with band B = D/L and thus we 
can directly express it as a linear combination of the functions gi.Ai^), with coefficients given directly by the samples 
of the product term at the points ri A (with n=:0,...,A^ — 1). This corresponds to what is done for example in 
Refjii, where in the final system of equations this term gives only a diagonal contribution consisting of its samples 
at the points n A. However, since we have assumed a band B for p{x), this assumption on the band of the product 
term cannot be verified (unless a constant potential energy function is considered) and thus, for high-order solutions, 
introduces a discrepancy with respect to the solutions from a Fourier analysis. 

In the following we discuss a better approach that that makes this direct-space method exactly equivalent to a 
Fourier method. 

We define M a generic odd (for numerical convenience) positive integer number (which accounts for the fact that 
in general it can be useful to exploit in the calculation more samples of the potential than those at the positions for 
which we want to evaluate the wave function p), and in the calculation we include the samples of the potential function 
f{x) at the N AI {= 2Q + 1, with Q a positive integer) points taken at intervals multiple of A' = L/{N M) = A/M. 
Assuming that the function f{x) is band limited with a maximum band B' = Q/L (in such a way that sampling it 
with an interval equal to A' the Nyquist criterion B' < 1/(2 A') is satisfied) and applying the sampling theorem, we 
can replace it with 
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Assuming for the wave funetion the expression ([51), we can therefore express the product term as: 
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(16) 
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In order to simplify the following calculations, we express also h{x) in terms of its samples. Since we have assumed the 
maximum bands B — D/L for p{x) and B' — Q/L for f{x), h{x) has at most a band equal to B" = B+B' ~ {D+Q)/L 
and thus for the sampling theorem we can write it in terms of its {M + 1)N — 1 samples in the period L, taken with a 
sampling interval A" = L/{{M + 1)N-1) = AN/{{M+1)N -1) which satisfies the Nyquist criterion B" < 1/(2 A"): 



The value of the samples is: 
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where we have defined as Si{n, m) and S2{n, £) the values of the g functions at the sampling points nA" . After some 
algebra, we obtain that 
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Finally, we need to consider only the projections of h{x) on the functions g^^A{x) (with fi ~ 0, . . . ,N ~ 1), i.e. to 
approximate h{x) with 
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In particular, we have that 
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where we have found that 
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In conclusion, substituting Eq. p^ into Eq. (|23p . we have found that 
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where /?£„ is the quantity inside the curl brackets, and the inner sum (with index m) inside the curly brackets depends 
only on n and can thus be computed once for every values of fx and £. Therefore the coefficients Pim can be easily 
computed with this expression and the formulas for 5*1, S2 and S3 that we have previously written. 



CORRESPONDENCE WITH THE FOURIER ANALYSIS 



Using this formulation, the method is exactly equivalent to the classical treatment of the problem in the transformed 
domain. 

In that case, we substitute the expressions of the wave function and of the potential function with their Fourier 
coefficients. In particular, we would compute the Fourier coefficients of the potential function by means of the DFT of 
M N samples of the function within the period L. In this way, we recast the problem into a system of linear equations 
having the Fourier coefficients of the wave function as unknowns. Then we reduce the (in principle infinite) matrix of 
the system to a finite size, disregarding the Fourier components above a fixed cut-off (spacial) frequency, both for the 
unknown wave function and for all the terms appearing in the equation. If we choose D/L as the frequency cut-off, 
we obtain an A^ x iV linear system (with N = 2D + 1). 

First of all, let us notice that if we write a periodic function z{x) in terms of its A^ = 2 Z? + 1 samples at the points 
n A within the period L 
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z{x)= Y<^A)gf,A{^) 



(26) 



we can compute its Fourier series coefficients Z^ noting that the Fourier coefficients [G£,a]p of the function gi^/^ix) 
are given by (exploiting the definition gi^A^x) as a repetition of sine functions) 
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Therefore the DFT of z{x) computed on TV = 2D + 1 samples corresponds to the exact Fourier series of the function 
with maximum band D / L that has the same samples as z{x) in the period L. 

In particular, this consideration can be applied to the potential function f{x): calculating its DFT on the M N 
samples (as we do in the Fourier analysis) corresponds exactly to substituting f{x) with the function f{x) of Eg. 1161 
(as we do in our approach). 

Moreover, limiting the Fourier components of the unknown wave function to frequencies that have a modulus less 
than D / L corresponds exactly to the assuinption that we have made working in the direct space that p{x) is band 
limited with band D/L, which has allowed us to use the sampling theorem with A^ = 2D + 1 samples. 

The same consideration is valid for the derivatives of p{x). The expressions that we have found for the derivatives 
can actually be found also starting from a Fourier analysis. Indeed, if p{x) is periodic and we assume it to be band 
limited with band B — D/L with Fourier coefficients i?p, we can write it and its derivatives in the form 
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The exponential function exp(i27rpa;/i) with \p\ < D can be seen as a band limited function with band B = D/L and 
thus expressed in terms of its N samples in the period L: 
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Substituting this expression inside Eq. (fSQ]) we obtain 
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Therefore (if we express Rp using Eq. (pS)) ) we have that 
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For s = 1 and s = 2 we have verified that this expression for ajg gives exactly the same results that are found working 
in the direct space (Eqs. (fTTj) and (jTSJ) ). 

Concerning the final treatment of the product term, in the reciprocal domain we consider only the frequency 
components of the product between —D/L and D/L, when operating a frequency cut-off. Since 1/(2A) = (2D + 
1)/(2L) = (D/L) + (1/(2L)) and all the terms are periodic with period L, and thus contain only frequency components 
multiple of l/L, this corresponds to limiting the band of the product term to the frequency range [— 1/(2A), 1/(2A)]. 
If H{f) is the Fourier transform of the function h{x) (/ being the spatial frequency), this corresponds to considering, 
in the Fourier domain, a spectrum Ho{f) — rect(A/)F(/), where rect(A/) is a function equal to 1 between — 1/(2A) 
and 1/(2A), and to outside this interval. The inverse Fourier transform of i?o(/) is equal to 
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Since the periodic function /io(a;) is band limited with band B < 1/(2A), we can express it in terms of its samples 
taken with sampling interval A: 
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(where we have exploited the fact that the sine function is even). This exactly corresponds with what we have done 
in the direct space. 

Therefore, with the approach we have described the solutions in the direct and in the reciprocal space are equivalent 
and correspond to linear systems of the same size. 

We notice that the periodic and band limited function p{x) can be equivalently expressed, using both the Fourier 
expansion and the sampling theorem, in the following two equivalent ways: 

D i2-Kp— '^^ 
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p=-D ^lU + i ^^^ 

where it is apparent that two different sets of orthonormal basis functions have been used: the functions 
exp(i27rpa;/L)/'\/2-D + 1 and the functions gi^^ix) (we have multiplied and divided by \/2D + 1 just to normal- 
ize the first set with respect to the scalar product @). Performing the solution in the reciprocal domain corresponds 
to using the first basis set, while performing the solution in the direct space corresponds to using the second basis 
set. The change from a basis to the other can be performed noting that each of the exponential functions appearing 
in the Fourier expansion (being band limited with band B < 1/(2A)) can be expressed in terms of its samples taken 
with sampling interval A in this way: 

: ~ > -Qj A(a;) = > Tirigi a{x) ■ (37) 

V2DTT ^ V2i? + 1 ' ^ "^ ' 

The matrix T containing Tip as {£,p) element allows to change the basis from the Fourier one to the direct space 
one. In particular, if the matrix of the direct-space linear system is M^ and the matrix of the reciprocal-space linear 
system is Mr, we have that Md = T Mr T~^ = T Mr T\ as we have numerically verified for the case described in the 
section IVH 

VI. APPLICATION TO THE SOLUTION OF THE DIRAC EQUATION IN AN ARMCHAIR 

GRAPHENE NANORIBBON 

Transport in graphene is described, within an envelope function approach, by the Dirac equation^. When consid- 
ering graphene nanoribbons, we have to enforce Dirichlet boundary conditions which couple the envelope functions 
corresponding to the two inequivalent Dirae points K and /\'. In particular, we take into consideration the solution 
of the Dirac equation in the case in which the potential inside the ribbon depends only on the transverse coordinate. 

As shown in RefsJ^ii^, the problem can be reformulated into a differential problem with periodic boundary condi- 
tions on a doubled domain. This new problem can be expressed in the following form: 

^di^ + ^^A(ll^ ^\W^ y\my) = -K^y) 

dy , (38) 

where the cr's are the Pauli matrices, (/?(?;) is a two-component function directly related to the four envelope functions 
of graphene and k^ is the longitudinal wave vector. Moreover, K is the modulus of the transversal coordinate of the 
Dirac points a,nd K = K - round {K/{tt/W)) tt/W is a quantity such that e'''^^^ = e-'"^^^ and that | A'| < n/W. 
Additionally, y is the transverse coordinate, W is the effective width of the ribbon (in the original problem, the 
Dirichlet boundary conditions had to be enforced at y = and y — W, where effective edges were located) and 
A(?;) = {U{y) — E)/^ is the potential term (U{y) is the potential energy, E is the Fermi energy, and 7 = hvp, where 
h is the reduced Planck constant and vp is the Fermi velocity in graphene). In this formulation, the problem has 
to be solved over the domain [0, 2 W] and the boundary condition corresponds to a periodic boundary condition on 
p{y) = 0{y)e-^''y. 
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FIG. 1. Potential energy considered in the armchair nanoribbon. 



Therefore we can rewrite the system in terms of the function p{y) as foUows (defining f{y) = K{W 

+ iKp{y) ] + cr^f{y)piy) = -KxPiy) 
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dy 
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The periodic boundary condition makes it possible to solve the problem in the Fourier domain or equivalently in 
the basis of g functions. Working in the basis of the g functions, we have rewritten the equation in the following way: 
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Projecting this equation (using the scalar product ^) onto the generic function gj^Aix) we obtain 
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These N ^ 2 D + 1 equations can be written in matrix form as an eigenproblem with eigenvalues —Kx and eigenvectors 
p{j A) which represent the values of the unknown function p{y) at the N points of the considered grid. 

In order to adopt the first approximate approach for the product term described at the beginning of the section HVl 
we have just to substitute /3ji with f{j A)Sji in Eq. (|42|) . 

In particular, wc have considered a potential U{y) (represented in Fig.[T]) equal to the sum of two Gaussian functions: 
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in a 8131 dimcr line wide (~ 1/^m wide) armchair nanoribbon, considering an electron injection energy E = 0.1 cV. 
In Figs. [2] and |3] we show the results that we have obtained considering A^ = 199 and M — 21 with the advanced 
approach in the direct space, and we compare them with those achieved by means of a Fourier analysis. In particular, 
in Fig. [5] we show the position of all the obtained longitudinal wave vectors k^ on the Gauss plane, and in Fig. [3] we 
represent the behavior of one of the four envelope functions corresponding to the eigenvalue with the largest real part. 
As we can see, the two methods yield exactly the same results. 

Instead, in Fig. 2] we show the values of the longitudinal wave vectors k^. obtained with the simplified treatment of 
the product term. As we can see from the highest order wave vectors, this method is not equivalent to the Fourier 
one, even though it gives quite good results for the low-order modes. 
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FIG. 2. Position on the Gauss plane of the all the longitudinal wave vectors Kx obtained with the Fourier analysis and with 
our advanced sine-based approach. 
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FIG. 3. Envelope function $a (y) corresponding to the eigenvalue with greater real part, computed with a Fourier analysis and 
with our advanced sine-based approach. 



VII. CONCLUSIONS 

We have proposed a method, based on a set of basis functions resulting from the periodic repetition of sine functions, 
for the numerical solution in the direct space of differential wave equations with periodic boundary conditions. This 
approach is fully equivalent (and requires the solution of a linear system with the same matrix size) to a Fourier 
analysis. We have applied this method to the solution of the Dirac equation in an armchair graphcne nanoribbon 
with longitudinal translational symmetry, verifying its equivalence to a reciprocal space solution. 
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